#################################################################
##DO SELF-REPORTING REGIMES MATTER? EVIDENCE FROM THE CAT########
## Cosette D. Creamer and Beth A. Simmons #######################
## February 2019 ###############################################
### Replication Code for In-Text Figures########################
### ########################
#################################################################
#################################################################
#################################################################

## Clear environment
rm(list=ls())
sink("CAT_In-Text Figures")
### Libraries
library(data.table)
library(MASS)
library(car)
library(SparseM)
library(ggplot2)

###### Figure 1: Trends in Protection Scores by Reporting Group###
#################################################################

df1 <- read.csv("figure1.csv", header=TRUE)

fig1 <- ggplot(df1, aes(x=Year, y=ValueFariss, group=Category)) + geom_line(aes(linetype=Category)) + theme_bw()
fig1 <- fig1 + labs(y="Average Human Rights Protection Score")
fig1
ggsave("Figure1.eps", width = 5, height = 3.5, units = "in")

#################################################################
#################################################################


################################################################
###### Figure 2: Estimated CmAT review treatment ###############
########## and treatment history effects ###### ################
#################################################################
########## Note: this uses output generated by ##################
########## running the "CAT_LHRPS MSM Analysis.R" script#########
#################################################################
###### Data loading and pre-processing #########################
## Results for LHRPS (t+2), Weighting Model 3####################
model2<-read.csv("CAT_LHRPS__model_2_alpha_0_iterations_20000.csv")


model2_point <- model2[1,]
model2_95ci_lower <- apply(model2, 2, function(x) quantile(x, .025))
model2_95ci_upper <- apply(model2, 2, function(x) quantile(x, .975))
model2_50ci_mean <- apply(model2, 2, function(x) quantile(x, .50))

out2low<-t(model2_95ci_lower)
out2upper<-t(model2_95ci_upper)
g2<-rbind(model2_50ci_mean,out2low,out2upper)
ci2<-t(g2)

results <- data.frame(method = c("Review at\ntime t", "Cumulative number of\nreviews at time t"),
                      point = c(model2_point$CATreview, model2_point$CATnumprevrevs),
                      ci95lower =c(model2_95ci_lower["CATreview"], model2_95ci_lower["CATnumprevrevs"]),
                      ci95upper = c(model2_95ci_upper["CATreview"], model2_95ci_upper["CATnumprevrevs"]))

######################## Plot point estimates
p = ggplot(results, aes(y=point, x=method))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=.5,linetype="dotted") 
p = p + geom_pointrange(aes(ymin=ci95lower,ymax=ci95upper),position="dodge", size= 1.5)
p = p + scale_y_continuous(name="Estimated average effect of one unit increase on Fariss score at time t+2", limits=c(-.6,.6))
p = p + scale_x_discrete(name="")
#p = p + geom_text(data=mentions_results_lags_pvals, aes(x=x, y=y, label=label), size=4, colour=c("black"))
# colour scheme
theme_bw1 <- function(base_size = 9, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + theme_bw1()
ggsave("Figure2.eps", width = 6, height = 2.5, units = "in")

#################################################################
###### Figure 3: Sensitivity Analysis for Review History Effects#
#################################################################
##############################################
alphas <- round(seq(-.17000, .170000, by=0.01),2)
pointestimates <- rep(NA, length(alphas))
upper_bounds <- rep(NA, length(alphas))
lower_bounds <- rep(NA, length(alphas))
upper_bounds2 <- rep(NA, length(alphas))
lower_bounds2 <- rep(NA, length(alphas))

for (i in 1:length(alphas)){
  alpha <- alphas[i]
  if (round(alpha) == alpha){
    alpha <- round(alpha)
  }
  simresults <- read.csv(paste("CAT_LHRPS__model_2_alpha_",alpha,"_iterations_20000.csv",sep=""))
  model2_point <- simresults[1,]
  model2_95ci_lower <- apply(simresults, 2, function(x) quantile(x, .025))
  model2_95ci_upper <- apply(simresults, 2, function(x) quantile(x, .975))
  
  model2_90ci_lower <- apply(simresults, 2, function(x) quantile(x, .05))
  model2_90ci_upper <- apply(simresults, 2, function(x) quantile(x, .95))
  
  
  pointestimates[i] <- model2_point["CATnumprevrevs"]
  upper_bounds[i] <- model2_95ci_upper["CATnumprevrevs"]
  lower_bounds[i] <- model2_95ci_lower["CATnumprevrevs"]
  upper_bounds2[i] <- model2_90ci_upper["CATnumprevrevs"]
  lower_bounds2[i] <- model2_90ci_lower["CATnumprevrevs"]
  
  print(alpha)
}



### Generate plot
postscript("Figure3.eps", width=4.5, height=4.5, horizontal = FALSE, 
           onefile = FALSE)
plot(x=alphas, y=pointestimates, type="l", lwd=3, xlim=c(-.17,.17), ylim=c(-.2, .5), xlab="Degree of Confounding: Y(Report) - Y(No Report)", ylab="Estimated marginal effect of one additional review", cex.lab=.7, cex.axis=0.7)
lines(x=alphas, y=upper_bounds, type="l", lty=2, lwd=2)
lines(x=alphas, y=lower_bounds, type="l", lty=2, lwd=2)
lines(x=alphas, y=upper_bounds2, type="l", lty=3, lwd=2)
lines(x=alphas, y=lower_bounds2, type="l", lty=3, lwd=2)
abline(v=0, lty =1,col="darkgrey")
abline(h=0, lty=3)
dev.off()


#################################################################
###### Figure 4: Shadow Reports ###### 
#################################################################
df2 <- read.csv("figure4.csv", header=TRUE)
df2$Report.Type <- factor(df2$Report.Type, levels = c("NHRI", "Domestic CSO", "International NGO"))
fig4 <-  ggplot(df2, aes(x=Year, y=value, fill=Report.Type))
fig4 <- fig4 + geom_bar(stat="identity") + theme_bw() + scale_fill_grey(name="Report Type")
fig4 <- fig4 + labs(y="Average shadow submissions, per report")
fig4
ggsave("Figure4.eps", width = 5, height = 3.5, units = "in")


#################################################################
######### Figure 5: Domestic Media Coverage of CAT ##############
#################################################################
df3 <- read.csv("figure5.csv", header=TRUE)
fig5 <- ggplot(df3, aes(x=Years, y=Value, fill=Reference))
fig5 <- fig5 + geom_bar(stat="identity", position=position_stack(reverse=TRUE)) + theme_bw() + scale_fill_grey()
fig5 <- fig5 + xlim(-6, 6) + labs(x = "Years Since Review", y= "Average Number of Articles")
fig5
ggsave("Figure5.eps", width = 5, height = 3.5, units = "in")

#################################################################
######### Figure 6: Legislative References in Chile ##############
#################################################################
df4 <- read.csv("figure6.csv", header=TRUE)
fig6 <- ggplot(data=df4, aes(x=Year, y=PercTotal, fill=Reference))+ geom_bar(stat="identity") + theme_bw() + scale_fill_grey(name="Reference") + ylab("% Total Sessions") + xlab("Year") + xlim(2004, 2018)
fig6 <- fig6 +  geom_vline(xintercept = 2004, linetype="longdash") + geom_text(df4, size=2, mapping=aes(x=2005, y=7.5, label="CmAT Review")) +  geom_vline(xintercept = 2009, linetype="longdash") + geom_text(df4, size=2,  mapping=aes(x=2011, y=7.5, label="CmAT Review"))  +  geom_vline(xintercept = 2007) + geom_text(df4, size=2, mapping=aes(x=2007, y=5.8, label="Report"))+  geom_vline(xintercept = 2017) + geom_text(df4, size=2, mapping=aes(x=2017, y=5, label="Report"))
fig6
ggsave("Figure6.eps", width = 5.5, height = 3.5, units = "in")


#################################################################
######### Figure 7: Legislative References in Argentina ##############
#################################################################
df5 <- read.csv("figure7.csv", header=TRUE)
fig7 <- ggplot(data=df5, aes(x=Year, y=PercTotal, fill=Reference)) + geom_bar(stat="identity") + theme_bw() + scale_fill_grey(name="Reference") + ylab("% Total Sessions") + xlab("Year")
fig7 <- fig7 +  geom_vline(xintercept = 2002) + geom_text(df5, size=2, mapping=aes(x=2001, y=12, label="Report")) +  geom_vline(xintercept = 2004, linetype="longdash") + geom_text(df5, size=2,  mapping=aes(x=2004, y=12, label="CmAT Review")) +  geom_vline(xintercept = 2015) + geom_text(df5, size=2, mapping=aes(x=2014, y=12, label="Report"))
fig7
ggsave("Figure7.eps", width = 5.5, height = 3.5, units = "in")

#################################################################
######### Figure 8: Legislative References in Mexico ##############
#################################################################
df6 <-read.csv("figure8.csv", header=TRUE)
fig8 <- ggplot(data=df6, aes(x=Year, y=PercTotal, fill=Reference)) + geom_bar(stat="identity") + theme_bw() + scale_fill_grey(name="Reference") + ylab("% Total Sessions") + xlab("Year")
fig8 <- fig8 +  geom_vline(xintercept = 2012, linetype="longdash") + geom_text(df6, size=2,  mapping=aes(x=2013, y=15, label="CmAT Review")) +  geom_vline(xintercept = 2011) + geom_text(df6, size=2, mapping=aes(x=2010, y=15, label="Report")) +  geom_vline(xintercept = 2006, linetype="longdash") + geom_text(df6, size=2,  mapping=aes(x=2006, y=15, label="CmAT Review")) +  geom_vline(xintercept = 2004) + geom_text(df6, size=2, mapping=aes(x=2003, y=15, label="Report"))
fig8
ggsave("Figure8.eps", width = 5.5, height = 3.5, units = "in")


sink(NULL)

